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1.  Foreword 

The  broad  mission  of  this  project  was  to  explore  new  classes  of  numerical  methods  that 
would  provide  new  approaches  for  modeling  complex  problems  in  elastostatics, 
elastodynamics,  and  wave  and  impact  problems.  We  believe  that  two  major  results 
were  obtained  in  the  course  of  this  work:  The  development  of  new  classes  of  so  called 
Discontinuous  Galerkin  Methods  (DGMs),  including  the  underlying  mathematical 
theory,  and  Generalized  Finite  Element  Methods  (GFEMs,  also  referred  to  as  hp-clouds, 
Partition-of-Unity  Methods  (PUMs))  for  linear  and  non-linear  boundary  value  problems. 
These  new  approaches  generalize  and  extend  so-called  mesh-free  techniques  and  are 
compatible  with  existing  finite  element  software.  At  the  same  time,  they  demonstrate 
significant  improvements  in  performance  over  traditional  finite  element  approaches. 
Summaries  of  work  done  under  each  of  these  topics  is  given  in  the  body  of  this  final 
report. 
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4.  Statement  of  the  Problem  Studied 

The  general  class  of  problems  studied  in  this  effort  were  boundary-and  initial-value 
problems  in  elastodynamics,  but  applications  to  broader  classes  were  also  considered. 
These  included  applications  in  convection-diffusion  phenomena,  gas  dynamics,  and 
viscous  incompressible  flow  as  modeled  by  the  Navier-Stokes  equations.  TTie  major  goal 
of  the  work  was  to  develop,  analyze,  and  implement  completely  new  types  of  methods 
which  could  overcome  numerous  shortcomings  of  existing  techniques  for  treating  large 
scale  simulations  in  solid  (and  fluid)  mechanics.  Two  avenues  of  research  were 
pursued,  one  involving  an  attempt  to  extend  our  earlier  work  on  Discontinuous 
Galerkin  Methods  for  hyperbolic  systems  to  problems  with  significant  physical  diffusion 
terms,  and  the  second  to  explore  new  methods  that  extend  and  make  competitive  ideas 
underlying  the  so-called  mesh-free  methodologies.  We  believe  that  efforts  in  both  of 
these  areas  were  extremely  successful.  Further  details  on  results  and  accomplishments 
are  given  in  the  following  paragraphs. 
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This  document  summarizes  results  obtained  on  a  project  aimed  at  developing  new  classes  of 
numerical  methods  for  the  analysis  of  problems  in  elastodynamics  and  elastostatics.  Two  significant 
classes  of  new  methods  were  developed,  analyzed,  and  implemented: 

1)  the  so-called  hp-Cloud  Method,  a  variant  of  the  meshfree  methods  built  on  partitions  of  unity 
generated  by  traditional  finite  elements  (also  referred  to  as  the  Generalized  Finite  Element  Method 
[GFEM])  and, 

2)  Discontinuous  Galerkin  Methods  for  broad  classes  of  transport  problems,  including  problems  with 
significant  diffusion.  These  new  methods  offer  numerous  advantages  over  traditional  schemes  for  a 
significant  class  of  applications.  A  summary  of  major  features  is  given  together  with  an  Appendix 
outlining  a  priori  error  estimates  and  convergence  proofs  for  various  Discontinuous  Galerkin  Methods. 
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5.  Summary  of  Most  Important  Results 

The  success  of  such  methods  as  finite  elements  on  a  broad  class  of  engineering  and 
scientific  applications  has,  in  many  ways,  made  obscure  several  serious  limitations  of 
these  approaches  in  treating  very  complex  problems  of  interest  in  applications  involving 
highly  nonlinear  phenomena  and  in  problems  which  require  high  fidelity  solutions.  Of 
these  limitations,  first  and  foremost,  is  the  a  long  standing  problem  of  developing 
computational  methods  for  transport  phenomena:  the  development  of  cell-wise 
conservative  schemes  that  can  deliver  very  high-order  local  accuracy.  The  search  for 
such  methods  has  been  "the  holy  grail"  of  computational  fluid  dynamics  for  over  two 
decades  and  has  lead  to  the  large  volume  of  papers  concerned  with  very  specialized 
techniques  for  hyperbolic  systems  that  employ  various  non-local  methods  to  get  high- 
order  accuracy.  Until  recently,  none  of  the  existing  techniques  had  sufficient  robustness 
or  generality  to  be  applicable  to  important  classes  of  three-dimensional  situations.  The 
new  DGMs  not  only  produce  local  cell-wise  conservative  approximative  schemes  but 
they  are  also  capable  of  delivering  high-order  polynomial  accuracy.  We  have  applied  the 
method  successfully  to  broad  classes  of  simulations  including  convection  diffusion 
problems,  three  dimensional  Euler-equations  for  compressible  flow,  2D  incompressible 
Navier-Stokes  equations,  and  we  are  in  the  process  of  implementing  the  algorithm  for 
large  scale  problems  in  impact  dynamics  and  penetration  mechanics. 

That  the  new  family  of  Discontinuous  Galerkin  Methods  is  applicable  to  problems  with 
diffusion  terms  was  announced  in  a  paper  by  Oden,  Bauman,  and  Babuska  [1]  published 
in  the  Journal  of  Computational  Physics.  This  particular  method  leads  to  an 
unsymmetrical  system  of  equations  for  the  diffusion  problem  which,  experimentally,  has 
proved  to  be  robust  and  has  the  ability  to  support  high-order  local  approximations. 
Indeed,  hp  versions  of  this  technique  have  been  developed  which,  with  adaptivity,  can 
yield  exponential  rates  of  convergence.  The  method  was  the  basis  of  a  dissertation  of 
Carlos  Bauman  and  his  version  of  the  method  has  become  known  in  the  literature  as  the 
Bauman-Oden  Method  (BOM).  As  noted  above,  the  BOM  has  subsequently  been 
extended  to  a  number  of  significant  two  and  three-dimensional  applications,  not  only  in 
computational  fluid  mechanics,  but  more  recently  in  treating  the  hyperbolic  coupled 
systems  encountered  in  the  imicible  displacement  equations  for  porous  medium. 
Extensions  of  the  BOM  to  these  problems  was  the  subject  of  the  dissertation  of  Beatrice 
Riviera  [2]  and  a  mathematical  analysis  of  some  of  its  properties  have  been  treated  in  a 
recent  publication  by  Riviera,  Gerault,  and  Wheeler  [3].  In  a  recent  paper,  summarized 
in  the  appendix,  proofs  of  convergence  and  a  priori  estimates  for  the  BOM  and  other 
DGMs  have  been  derived.  We  also  review  in  the  appendix  to  this  report  details  of  the 
formulation  and  a  comparison  of  this  method  with  other  interior  penalty  methods  and 
discontinuous  methods  that  have  been  subsequently  proposed.  Special  features  of  the 
method  are  listed  as  follows: 

1.  Local  element  polynomial  approximations  of  arbitrary  order  are  constructed  over 
each  cell  providing  a  basis  for  local  spectral  type  approximations. 

2.  Continuity  conditions  on  the  solution  values  and  the  fluxes  at  cell  boundaries  are 
enforced  in  a  weak  form. 

3.  A  posteriori  error  indicators  have  been  derived  and  used  in  an  adaptive  procedure 
that  has  allowed  the  implementation  of  the  method  on  non-uniform  hp  meshes;  a 
number  test  problems  have  been  treated  which  demonstrate  exponential 
convergence. 

4.  Most  importantly,  the  scheme  is  cell-wise  conservative,  meaning  that  fluxes  are 
balanced  on  each  element;  as  far  as  is  known  by  the  investigators,  this  represents  the 
first  high-order  (greater  than  first  or  second  order)  conservative  FEM  scheme. 


5.  A  priori  and  a  posteriori  error  estimates  have  been  derived  for  linear  applications 
including  second  order  elliptic  problems,  convection  diffusion  problems  in  two 
dimensions,  and  first  order  hyperbolic  problems  in  two  dimensions.  The  estimated 
convergence  rates  have  been  confirmed  in  several  numerical  experiments.  Tests 
indicate  for  p-versions  of  the  method,  sub-optimal  rates  are  achieved. '  Since 
publication  of  our  work,  a  number  of  new  variations  of  the  method  have  been 
proposed  by  other  authors  which  include  the  use  of  penalized  or  stabilizing  terms 
which  may  lead  to  robustness  of  the  scheme  and  improve  rates  of  convergence,  but 
may  also  destroy  the  conservation  properties  of  the  method. 

6.  Three  dimensional  versions  of  the  BOM  have  been  developed  for  applications  in 
computational  fluid  dynamics,  (as  noted  above,  in  particular,  Euler-equations  and,  in 
two  dimensions,  the  incompressible  Navier-Stokes  equations).  These  have  been 
applied  to  a  significant  number  of  benchmark  problems. 


We  believe  these  methods  may  have  important  implications  for  coupled  problems  in 
nonlinear  dynamics.  Work  remains  to  be  done  in  deriving  sharp  a  posteriori  error 
estimates  of  these  nonlinear  problems,  including  local  estimates  for  quantities  of  interest. 


Mesh-free,  Partition  of  Unity,  and  Generalized  Finite  Element  Methods 

Another  area  explored  during  the  course  of  this  project  was  the  development  of  new 
families  of  schemes  for  treating  both  linear  and  nonlinear  boundary  value  problems  in 
mechanics  and  physics.  Our  work  began  on  this  subject  in  the  mid-1990's  with  the 
development  of  so-called  hp  Clouds.  These  were  techniques  based  on  the  use  of 
Moving-Least-Squares,  (MLS)  techniques  for  deriving  basis  these  functions  for  finite 
dimensional  spaces,  functions  are  built  on  overlapping  domains  which  supported 
polynomials  of  arbitrary  degree.  In  our  earlier  work,  we  were  able  to  develop  h,  p,  and 
hp  versions  of  these  techniques  and  to  develop  adaptive  versions  of  these  schemes  that 
demonstrated  exponential  convergence. 

Subsequent  work  showed  that  while  these  schemes  could  exhibit  extraordinarily  high 
rates  of  convergence  they  were  not  competitive  with  existing  finite  element  and  finite 
volume  techniques  because  of  the  significant  expense  required  in  quadrature  of  the  basis 
functions  defined  on  spherical  domains.  This  proved  to  be  expensive  and  cumbersome 
and  required  a  computational  effort  that  dominated  the  entire  cost  of  their 
implementation.  In  general,  we  have  concluded  that  the  cost  of  generating  a  partition  of 
unity  method  on  the  irregular  domains  using  the  methods  of  least  square,  a  common 
step  in  most  of  the  mesh-free  methods  now  in  use,  is  prohibitive. 

In  1997  and  1998,  we  began  exploring  the  use  of  conventional  finite  element  methods  as 
a  device  to  generate  a  partition  of  unity  over  a  given  computational  domain.  Such  an 
approach  automatically  overcame  the  quadrature  problem,  since  quadrature  points 
were  already  embedded  in  the  finite  element  master  elements,  but  they  (FEMs)  provided 
nodal  basis  functions  on  overlapping  domains  on  which  high-order  polynomial 
approximations  and  other  types  of  Treffz  approximations  could  be  constructed.  While 
the  resulting  methods  could  not  qualify  as  "mesh-free,"  they  possess  a  number  of  very 
desirable  properties  which  ultimately  proved  to  be  superior  to  conventional  finite 
element  techniques.  We  have  termed  these  methods  GFEMs:  Generalized  Finite 
Element  Methods.  They  were  first  reported  in  papers  by  Oden,  Duarte,  and 
Zienkiewicz,  and  in  a  one-diniensional  case,  by  Babuska  and  Melenik  [4].  With  the 
addition  of  special  functions,  these  techniques  can  provide  accelerated  rates  of 


convergence  in  cases  where  traditional  finite  element  methods  diverge  or  converge  very 
slowly.  These  include  problems  in  which  one  has  very  rough  highly  oscillatory 
coefficients,  boundary  layers,  interfaces,  singularities,  and  other  local  features  that 
influence  convergence. 

One  of  the  most  important  applications  of  the  GFEM  has  been  in  the  treatment  of  crack 
problems  and  propagation  of  cracks  through  domains.  We  have  developed  techniques 
for  allowing  cracks  to  occur  in  a  given  domain  and  to  propagate  through  the  domain 
without  altering  the  original  mesh.  Commercial  engineering  software  companies  have 
expressed  interest  in  our  results  and  one  (Altair  Engineering)  has  used  these 
methodologies  in  applications  involving  two  and  three-dimensional  crack  propagation 
problems.  Apparently,  our  techniques  have  formed  the  basis  of  certain  versions  of 
commercial  products  that  are  or  soon  will  be  on  the  market. 

In  many  of  the  applications  and  experimental  tests  of  these  GFEMs,  performances  have 
been  observed  which  are  considerably  faster  and  more  efficient  than  traditional  finite 
element  methods.  As  noted  earlier,  one  can  also  cite  problems  in  which  these  methods 
exhibit  extraordinarily  high  convergence  rates  while  FEMs  do  not  converge  at  all.  The 
subject  of  GFEM  has  been  further  explored  in  by  dissertation  by  Kevin  Copps  [5];  it  has 
been  the  subject  of  a  NSF  workshop  and  is  a  popular  topic  in  conferences  on 
computational  science  and  engineering.  The  full  exploitation  of  these  methods  in 
treating  significant  non-linear  problems  remains  to  be  explored.  We  feel  that  one  of  the 
most  important  features  of  GFEMs  over  some  of  the  recently  reported  mesh-free 
methods  is  that  they  can  be  easily  embedded  into  existing  finite  element  codes.  The 
GFEM,  in  other  words,  can  be  used  to  add  to  existing  finite  element  packages  the  ability 
to  treat  a  long  list  of  special  features,  particularly  the  propagation  of  cracks  across  fixed 
meshes  and  the  treatment  of  singularities  in  boundary  layer  effects.  We  have  derived  a 
priori  and  a  posteriori  error  estimates  for  these  methods  and  have  developed  adaptive 
versions  of  them. 

One  area  relative  to  GFEMs  that  remains  to  be  more  thoroughly  explored  is  the 
development  of  preconditioners  for  the  method  for  very  large  problems,  including  the 
use  of  the  techniques  in  a  parallel  computer  environment.  Since  the  GFEM  technique 
involves  superimposing  polynomial  type  approximations  (or,  in  some  cases,  special 
functions)  on  a  model-based  partition  of  unity  generated  by  traditional  FEMs,  in  most 
cases  rank  deficient  stiffness  matrices  are  produced.  Thus,  special  techniques  must  be 
developed  to  solve  the  resulting  linear  systems.  We  have  developed  and  implemented 
one  technique  that  seems  to  work  satisfactorily,  but  a  great  deal  of  additional  theoretical 
work  needs  to  be  done  to  develop  fully  robust  iterative  techniques  for  handling  the 
linear  systems  generated  by  GFEMs.  We  believe  this  problem  is  solvable,  and  certainly 
is  an  area  worth  further  study. 
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APPENDIX 


1.  Introduction 

There  has  been  renewed  interest  in  Discontinuous  Galerkin  Methods  (DGM)  re¬ 
cently,  primarily  due  to  the  discovery  that  variants  of  these  methods  could  be  used 
effectively  to  solve  diffusion  problems  as  well  as  problems  of  pure  convection.  One 
such  DGM  was  presented  in  the  dissertation  of  Baumann  [7]  and  reported  in  the 
paper  of  Oden,  Babuska,  and  Baumann  [18];  summary  of  other  versions  of  DGMs 
and  a  lengthy  historical  review  of  this  subject  can  be  found  in  the  record  volume 
edited  by  Cockburn,  Karniadakis,  and  Shu  [10].  The  DGM  possesses  a  number 
of  important  properties  that  set  them  apart  from  traditional  conforming  Galerkin- 
finite  element  methods:  they  are  elementwise  conservative,  can  support  high  order 
local  approximations  that  can  vary  nonuniformly  over  the  mesh,  are  readily  paral- 
lelizable,  and,  for  time-dependent  problems,  lead  to  block-diagonal  mass  matrices, 
even  for  high-order  polynomial  approximations.  These  properties  make  DGMs 
attractive  candidates  for  a  broad  collection  of  applications. 

Several  papers  have  been  published  in  the  mathematical  literature  on  a  priori  error 
estimates  for  various  DGMs  for  diffusion  problems.  In  particular,  an  analysis  of 
one-dimensional  versions  of  the  Baumann-Oden  method  was  reported  by  Babu§ka, 
Oden,  and  Baumann  [2].  Error  estimates  for  several  types  of  DGMs  and  for  the 
related  Internal  Penalty  Galerkin  Methods  were  presented  in  the  dissertation  of 
Riviere  [19]  and  in  the  paper  of  Riviere,  Wheeler,  and  Girault  [20].  Several  other 
studies  on  a  priori  error  estimates  for  DGMs  have  appeared  recently;  see,  for  exam¬ 
ple,  the  report  of  Chen  [9]  and  the  analysis  of  Siili,  Schwab,  and  Houston  [22,15]. 
Convergence  analysis  of  other  variants  of  DGM  can  be  found  in  [10]. 

Here  we  present  a  detailed  derivation  of  a  priori  error  estimates  for  several  ftp- 
versions  of  DG-finite  element  methods  for  linear  diffusion  problems  (the  Poisson 
problem)  on  two-dimensional  domains.  In  some  cases,  important  steps  in  our  anal¬ 
ysis  follows  the  approach  of  Riviere,  Wheeler,  and  Girault  [20],  but  other  steps 
differ  in  detail.  We  present  a  series  of  approaches  in  which  different  versions  of 
DGMs,  including  those  with  penalty  terms,  can  be  analyzed.  Our  final  estimates 
differ  in  predicted  rates  of  convergence  with  respect  to  the  polynomial  degree  p 
obtained  in  [20,19]  and  reflect  rates  consistent  with  the  computed  results  of  Bau¬ 
mann  [7]. 
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2.  Notations  and  Preliminaries 

We  shall  choose  the  domain  as  a  bounded  open  set  in  with  Lipschitz  con¬ 
tinuous  boundary  We  will  denote  r©  the  part  of  die  boimdary  dii  on  which 
Dirichlet  conditions  are  prescribed  and  the  part  on  which  Neumann  conditions 

are  prescribed.  Formally,  the  boundary  9Q  is  decomposed  into  the  parts  To  and 
such  that  Fp  U  =  dQ,  and  Fp  D  F^  =  0. 


2.1.  Finite  Element  Partition 

Let  !Pf,  denote  a  partition  of  the  domain  Q,  i.e.  2^  is  a  finite  collection  of  Ng  open 
subdomains  (elements)  K,-,  i  =  1, 2, . . . ,  Ng,  such  that: 

U  =  i#/. 

K,e25, 

The  size  and  shape  of  an  element  K,,  or  simply  K,  of  2^  are  measxured  in  terms  of 
two  quantities,  and  ptc,  defined  as: 

hx  =  diam(K), 

Pk  =  sup  {diam(®);  ®  is  a  ball  contained  in  K}. 

We  also  introduce  tiie  parameter  h  associated  with  the  partition  2^: 


h  =  maxhK-  (2-1) 

Ke2}, 

Definition  A  family  (2),}  of  partitions  %  is  said  to  be  shape  regular  as  h  tends  to  zero  if 
there  exists  a  number  q>Q,  independent  ofh  and  K  such  that: 

^<Q,  'iKen.  (2.2) 

PK 

In  this  appendix,  all  partitions  25,  are  assumed  to  be  shape-regular. 

In  addition,  we  shall  associate  with  each  element  K  the  element  boundary  dK.  The 
unit  normal  vector  outward  from  K  (resp.  K,)  is  denoted  by  n  (resp.  n|,). 

Given  a  partition  2},,  we  shall  denote  the  collection  of  edges  of  2},  (points  in  one 
dimension,  faces  in  three  dimensions)  by  the  set  =  {7/},  Z  =  1, . . .  ,  N.y.  Edges 
represent  here  open  subsets  of  either  Q  or  dFl.  We  thus  introduce  the  set  of 
interior  edges  as: 

If«f=  (2-3) 
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so  that: 

N7  _  _  _ 

U  7/  =  Td  u  r^/  u  Pint. 

1=1 

In  the  same  way,  we  shall  decompose  'Ef,  into  three  subsets  as: 

~  '^h,D  G  ‘^h,N  G  ^h,int- 

Then,  7  €  if  it  lies  on  To,  and  7  €  if  it  lies  on  r^.  Moreover,  as  shown 
in  Fig.  1,  '■fij  6  E^^nt  denotes  an  edge  (interface)  between  two  adjacent  elements  K,- 
and  Kj,  where  by  convention  i  >  j.  For  each  edge  7,  we  also  associate  a  unit  normal 
vector  n.  In  the  case  7  is  an  edge  associated  with  an  element  K,  adjacent  to  dQ,  i.e. 
7  €  U  ,  the  unit  normal  vector  is  simply  defined  as  n  =  n|,-.  For  an  interior 
edge  jij  G  with  the  convention  i  >  j,  n  is  chosen  as  the  unit  normal  vector 

outward  from  Ki,  so  that  n  =  n|,  =  -n\j  (see  Fig.  1).  In  subsequent  analyses,  C  will 
denote  generic  positive  constants,  not  necessarily  the  same  in  different  places. 

Remark  1  Using  simple  geometrical  properties,  one  can  show  that  each  edge  7  in  a  shape- 
regular  partition  satisfies: 

-hK  <Pk<  I7I  <  (2.4) 

Q 

where  |7|  denotes  the  length  of  In  other  words,  hf^  and  7  are  equal  within  a  constant. 
Therefore,  we  will  interchangeably  use  hic  or  7  (preferably  hfc). 


2.2.  Spaces 

Let  s  be  a  positive  integer.  For  any  given  open  set  S  (S  may  define  the  whole 
domain  Q,  an  element  K  of  E;,,  or  an  edge  7  of  2),),  the  spaces  H®(S)  will  denote  the 
usual  Sobolev  spaces  with  norm  ||  •  ||.s,s.  hi  the  particular  case  in  which  S  represents 
Q,  the  norm  will  simply  be  denoted  IHIs.  Moreover,  H^(S)  is  the  set  of  functions  in 
H^(S)  which  vanish  on  the  boundary  dS  of  S,  i.e. 

Ho(S)  =  {v€  H^(S);v  =  0  on  5S}, 

and  H(div,  S)  denotes  the  space: 

H(div,  S)  =  {v  G  (L2(S))2;  V  •  v  G  l2(S)}. 

The  so-called  (mesh-dependent)  broken  space  will  be  defined  as: 

W(Eh)  =  {ve  L2(Q);  i;|k  G  H^(fC),  VfC  G  E^}. 
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Figure  1.  Element  interface  7,)  and  unit  nor¬ 
mal  vector  n. 

The  norm  associated  with  the  space  is  given  as: 

ll-lkA  =  f  S  ll'’ll?,«') 

\»:e«  / 

where  ||y||s,jR:  is  the  Sobolev  norm  on  K. 

We  will  consider  finite  element  spaces  of  polynomial  functions,  possibly  dis¬ 
continuous  at  the  element  interfaces,  such  as: 

=  {ve  L2(Q);  o  f)  e  Ppih  VfC  6  Tj,}  (2.5) 

where  Fr  is  the  affine  mapping  from  the  master  element  K  to  the  element  K  in  the 
partition,  and  Pp(K)  is  the  space  of  polynomial  functions  of  degree  at  most  p  on  K. 

In  hp  methods,  the  polynomial  degree  can  actually  vary  from  one  element  to  the 
other.  Denoting  pK  the  pol5momial  degree  associated  with  the  element  fC,  we  define 
the  global  value  p  for  the  partition  fPfi  as: 

p  —  min  pk-  (2.6) 

One  advantage  of  DGMs  over  conventional  hp  finite  element  methods  is  that  the 
polynomial  degrees  pjc  do  not  necessarily  match  at  the  interfaces  of  the  elements. 
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3.  Formulations  for  the  Poisson  Model  Problem 


3.1.  Model  Problem 

We  shall  consider  here  the  following  Poisson  model  problem:  find  the  scalar  func¬ 
tion  u  which  is  the  solution  of 

—Am  -I-  cm  =  /,  in  Q,  (3-1) 

and  which  satisfies  the  boundary  conditions: 


M  =  Mo,  onFo, 
n-Vu  =  g,  onF^. 


(3.2) 


Here  /  e  L^{Q)  represents  the  load  scalar  and  c  is  a  positive  constant  over  the  do¬ 
main  Q. 

We  now  proceed  with  ttie  derivation  of  weak  formulations  of  the  Poisson  equation 
on  discontinuous  spaces.  Let  u,  for  the  moment,  be  a  sufficiently  smooth  func¬ 
tion.  The  regularity  of  m  shall  be  discussed  later  in  the  appendix,  namely  in  Sub¬ 
section  3.3.  Multiplying  (3.1)  by  a  function  v  in  H\Th)  and  integrating  over  the 
domain  Q,  we  obtain: 


•  Vu  +  cu)  vdx  = 


Unlike  the  classical  continuous  finite  element  approach,  we  shall  first  decompose 
the  integrals  in  the  above  equation  into  element  contributions 


y  -  f  (V-VM)udx-h  y  f  cuvdx=  X  / 
and  then  integrate  by  parts,  so  that: 

y  f  (Vu-Vv  +  cuv)dx-  y  f  {n-Vu)vds=  X  f  (3.3) 

^  K62), 


We  observe  that  the  boundary  integrals  are  defined  on  each  element  boundary; 
those  are  now  splitted  according  to  the  type  of  boundary  such  as: 


X  f  {r\-'Vu)vds  —  X  f  (n-VM)uds 

X  /  (n.VM)i;ds 

-p  X  /  (n- VM),  t;,  -l-(n.  VM)yt;^ds. 
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where  Vi  and  Vj  denote  tire  restrictions  of  v  on  the  elements  Kj  and  Kj  respectively. 
In  the  same  way,  (n  •  Vm),-  and  (n  •  Vu)j  represents  the  restrictions  of  the  flux  n  •  Vw 
on  Ki  and  Kj. 

In  general,  except  occasionally  to  avoid  confusion,  we  shall  simplify  the  notation 
of  dxese  boimdary  integrals,  by  rewriting  them,  for  instance. 


2)  f  {n-Vu)vds  = 
^  f  (n-'Vu)vds  — 


Moreover,  the  treatment  of  the  interior  boundary  integrals  is  as  follows.  Given  an 
edge  jij  €  £h,M  shared  by  two  adjacent  elements  fC,-  and  Kj,  i  >  j,  we  first  note  that: 


(n  ■  Vm),-  Vi  +  (n  •  Vm)^  Vj  =  n-  (Vu)i  Vj-n-  (yu)j  Vj, 

where  n  is  now  the  unit  normal  vector  with  respect  to  the  edge  7,^  as  defined  in  the 
previous  section.  By  analogy  with  the  formula  below  where  a,  b,c  and  d  are  real 
numbers: 

ac-bd  =  |(fl  +  b)ic  -d)  +  |(a  -  b){c  +  d),  (3.4) 

we  can  write  the  integrand  as: 


n  •  (Vm),-  Vi -ti-  (Vu)j  Vj 

=  i  (n  •  (Vm),-  +  n  •  (Vm)^)  (u,-  -  Vj)  +  ^  (n  •  (Vm),-  -  n  •  (Vm)J  (t;,-  +  Vj) 
=  (n  •  Vm)  [v]  +  [n  •  Vm]  {v)  . 


Here  [v]  and  {v)  respectively  denote  the  jump  and  average  of  on  an  interior  edge 
lijf  i  >  //  of  function  v  G  H^(K{)  x  H^{Kj),  s  >  1/2,  i.e. 

[z;]  =  Vi  -Vj, 

1 

(v)  =  -iVi  +  Vj). 

We  conveniently  extend  the  definition  of  [zj]  and  (v),  following  Chen  [9],  to  an  edge 
7  lying  on  Fd  as: 

[Z7]  =  V, 

(v)  =  V. 
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It  allows  us  to  combine  the  interior  and  Dirichlet  boundary  terms  in  only  one  inte¬ 
gral  as: 


^  f  {ii  ■  Vu)i  Vi +  {n-Vu)jVjds+  ^  j  {n-Vu)vds 


=  f  (n- VM)[z;]-f  [n- Vu](i;)ds. 


Remark  2  Note  that  when  u  e  the  fluxes  [n  •  Vu]  are  continuous  almost  every¬ 

where  in  Q,  which  yields 

f  [n  •  Vm]  (u)  ds  =  0,  VveH\lPh).  (3.5) 


Consequently,  (3.3)  can  now  be  reduced,  when  u  €  H^(Q.)  and  applying  the  Neu¬ 
mann  boundary  condition,  to: 

y,  f  {Vu-'Vv-\-cuv)dx-  f  (n-VM)[u]ds=  Y  f  fvdx+  [  gvds. 

Kelj,  ■'hniuro  KeiPi, 


We  introduce  the  following  bilinear  form  B(-,  •)  defined  on  x  H^{(Ph)  and  the 

linear  form  L(-)  defined  on  H\lPh)  such  as: 

B(u,v)=  2]  /  (Vm  •  Vu-!-CMt;)dx,  (3.6) 

K625, 

f(t^)=  X  [  fvdx+  f  gvds.  (3.7) 

We  also  consider  the  bilinear  form  /(•,  •)  on  x  H^{lPh),  which  incorporates  all 

boundary  integrals  on  Fjnf  and  F p,  as: 

J{u,v)=  f  (n-Vu)[u]ds.  (3.8) 

Then,  a  general  discontinuous  weak  formulation  of  the  Poisson  equation  reads: 

Biu,v)  -  J{u,  V)  =  F(u),  Vu  e  H^lPh).  (3.9) 

This  above  variational  form  constitutes  the  starting  point  to  derive  formulations  of 
various  Discontinuous  Galerkin  Finite  Element  Methods  (concisely,  DGMs.) 
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3.2.  Weak  Formulations  and  Finite  Element  Discretizations 

All  the  formulations  presented  below  use  the  observation  that,  for  u  €  H^(C2)  n 
H\Th),  the  jump  [«]  vanishes  on  each  7^-: 

f  u[M]ds  =  0,  'ivGL^'yij).  (3.10) 

Jyij 

It  follows  that: 


f  (n  •  Vv)  [u]  ds  =  0,  Vi;  e  H^^Ph). 

•'^ni 


(3.11) 


Moreover,  the  Dirichlet  boundary  condition  can  be  applied  in  the  following  weak 
manner: 

f  (n:Vv)uds=  f  (n  •  Vv)  uq  ds,  Vu  G  H\‘Ph).  (3.12) 

JTd  JTd 

Therefore,  introducing  the  linear  form  /o(-)  defined  as: 


j^(v)  =  /  (n  •  Vv)  «o ds,  Vt;  G  (3-13) 

JTo 

we  observe  that,  for  u  G  H^(Q)  D  ®*^d  m  =  Mo  on  Fq, 

}{v,  u)  =  Joiv),  Vu  G  (3.14) 


3.2.1.  Global  Element  Method  -  GEM 

Introducing  the  bilinear  form  ®_(-,  •)/  the  subscript  -  referring  to  the  fact  that  we 
substract  the  term  J(v,  u)  to  the  left  hand  side  of  (3.9),  and  the  linear  form 

'S-{u,v)  =  B(u,v)-Jiu,v)-J(v,u), 

!F-(v)  =  Fiv)-h{v), 

the  Global  Element  Method  consists  in  finding  u  such  that: 

(M,  z;)  =  fF-{v),  Vi;  G  H^Ph)-  (3-16) 

One  advantage  of  this  method  is  that  it  defines  a  symmetric  problem.  On  the  other 
hand,  a  sigiuficant  disadvantage  is  that  the  bilinear  form  is  not  guaranteed  to  be 
semi-positive  definite.  When  dealing  with  time-dependent  problems,  this  could 
imply  that  some  eigenvalues  have  negative  real  parts,  causing  the  formulation  to 
be  unconditionally  unstable. 

The  corresponding  finite  element  discretization  of  the  above  problem  consists  in 
finding  G  such  that; 

'B-{Uh,v)  =  f-iv),  'iv  G  (3.17) 

This  method  was  introduced  by  Delves  et  al  [11-14]  with  the  particular  objective 
of  accelerating  convergence  of  iterative  schemes. 
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3.2.2.  Symmetric  Interior  Penalty  Galerkin  Method  -  SIPG 

To  enforce  stability  of  the  discontinuous  method,  i.e.  continuity  of  the  solution  at 
the  interface  of  the  elements,  penalty  terms  have  been  added  to  the  formulation  by 
Arnold  [1]  and  Wheeler  [23].  Let  us  introduce  the  following  penalty  terms: 

}^(u,v)=  y  /  a[u][v]ds+  y  auvds=  /  a[u][v]ds, 

and 

y  /  (TUovds=  /  (TUovds, 

7eE»,D‘'T' 

where  cr  represents  the  penalty  parameter  which  depends  on  the  length  of  the 
edges  jij  and  7  and  the  polynomial  degree  used  in  the  elements;  namely  a  = 
a{h,  p).  Then  the  SIPG  method  is  similar  to  the  GEM  except  for  the  penalty  terms. 
Indeed,  introducing  the  forms: 

^(u,v)  =  B(u,v)-J(u,v)-J(v,u)-l-J^(u,v), 

!Ff{v)=Fiv)-jo(v)  +  JS{v), 

the  Symmetric  Interior  Penalty  Galerkin  problem  is  to  find  u  such  that: 

^{u,v)  =  ^veH\Th).  (3.19) 

Note  that  when  a  takes  on  the  value  zero,  we  naturally  retrieve  the  GE  method. 
The  finite  element  analogue  of  problem  (3.19)  is  to  find  «/,  G  such  that: 

‘BZiUh,  V)  =  (0),  \fv  G  (3.20) 

Remark  3  Following  Baker  and  Karakashian  [5,6,16],  we  consider  a  variant  of  the  SIPG 
method.  Instead  of  using  the  formula  (3.4),  one  may  use: 

ac  —  bd  =  ac-ad  +  ad  —  bd  =  a{c  —  d)  +  (a  -  b)d  (3.21) 

so  that,  by  analogy: 


n  •  (Vh),  Vj-n-  (yu)j  uy  =  n •  (Vw),- [u]  +  [n  •  V«] Vj 


and,  since  the  fluxes,  for  u  G  are  continuous  across  the  interelement  boundaries,  we 

have: 
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The  new  bilinear  form  for  the  boundary  terms  is  now  defined  as: 

Iiu,v)=  f  n-iVu)i[v]ds 

JTintUTD 

so  that  the  new  formulation  reads:  Find  u  €  HH^)  n  H\Th)  such  that,  for  all  v  €  Ffi(tPh), 

B{u,v)-Iiu,v)-I{v,u)  +  r{u,v)  =  F{v)-Jo{v)  +  Joiv).  (3.22) 

We  now  see  that  we  recover  the  SIPG  method  from  the  Baker-Karakashian  formulation  by 
replacing  the  term  n  •  (Vm),'  by  (n  •  Vm).  It  follows  that  all  the  properties  associated  with 
the  SIPG  method  will  also  apply  to  the  Baker-Karakashian  formulation. 


3.2.3.  Discontinuous  hp  Galerkin  FE  Method  -  DGM 

The  discontinuous  Galerkin  method  by  Baumann  et  al.  [7,18]  differs  from  the  Global 
Element  Method  by  just  a  sign.  Indeed,  by  introducing  the  forms: 

‘B+(u,v)=B(u,v)-Ku,v)  +  Kv,u), 

:r+iv)  =  F{v)-Fjo{v), 

the  DG  formulation  reads:  Find  u  such  that 

®+(M,  t;)  =  f+iv),  e  H\lPf).  (3.24) 

It  is  straightforward  to  show  that  ihe  bilinear  form  is  positive  semidefinite. 

The  associated  finite  element  version  of  the  DG  method  consists  then  in  finding 
such  that 

‘B+iUh^v)  =  %.{v)^  Vu  G  (3.25) 


3.2.4.  Non-Symmetiic  Interior  Penalty  Galerkin  Method  -  NIPG 

This  method  was  introduced  by  Riviere  [19]  and  Suli,  Schwab  and  Houston  [22,15] 
and  is  inspired  from  the  DG  method  with  the  addition  of  penalty  terms.  The  new 
bilinear  and  linear  forms  read: 

Bfiu, v)  =  B(u,  v)  -  Ku,v)  +  J{v,  u)  +  r{u,v), 

n(P)  =  F{v)  +  h{v)  +  lQ{'of 

so  that  the  problem  to  solve  by  the  NIPG  method  becomes:  Find  u  such  that 
Bfiu,v)= 

Once  again,  we  may  consider  DG  as  a  special  case  of  NIPG  with  <t  =  0. 


(3.27) 
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The  finite  element  problem  corresponding  to  the  NIPG  formulation  (3.27)  is  to  find 
Uh  €  such  that 

,  V)  =  \/v  e  (3.28) 

The  four  methods  presented  thus  far  are  all  very  similar,  except  for  a  plus  or  mi¬ 
nus  sign  in  front  of  the  term  J{v,  u)  and  the  addition  of  a  penalty  term  v)  or 
not.  We  shall  now  see  how  these  changes  modify  the  properties  of  the  respective 
formulations. 


3.3.  Equivalence  of  Strong  and  Weak  Problems 

We  shall  show  the  equivalence  of  the  strong  and  weak  formulations  only  with  re¬ 
spect  to  the  Global  Element  method.  The  results  are  identical  for  the  other  formu¬ 
lations,  namely  the  SIPG,  DG  and  NIPG  methods.  Existence  of  solutions  of  the 
discontinuous  formulations  is  then  somewhat  guaranteed.  However,  we  empha¬ 
size  here  that  Theorem  3.1  does  not  infer  anything  about  the  uniqueness  of  the 
solutions.  This  question  still  remains  an  open  issue. 

Theorem  3.1  (GE  Method)  Let  u  e  C^(Q)  be  the  solution  of  Problem  (3.1)-(3.2).  Then 
u  satisfies  the  weak  formulation  (3.16).  Conversely,  if  u  E  H^(Q)  D  if'ifPf)  is  a  solu¬ 
tion  of  (3.16)  then  u  satisfies  the  partial  differential  equation  (3.1)  and  boundary  condi¬ 
tions  (3.2). 


Proof:  The  first  part  of  the  theorem  has  been  proved  along  with  th^erivation  of 
the  Global  Element  formulation,  since  (3.9)  is  satisfied  when  u  €  C^(ii). 

The  converse  follows  the  proof  given  in  Riviere  [19].  Let  (D(K)  C  H^{K)  be  the 
space  of  infinitely  differentiable  functions  with  compact  support  on  element  K  and 
let  V  e  ID(K).  Then  (3.16)  gives: 


/  {Vu  ■Vv  +  cuv)dx=  /  fvdx 
Jk  Jk 


which  implies,  after  integration  by  parts  and  since  v  is  arbitrary  in  ^I>(K),  that 


—Am  +  cu  =  f,  a.e.  in  K. 


(3.29) 


Next,  we  consider  an  interior  edge  jij  shared  by  the  elements  Ki  and  Kj.  Let  v  be 
a  function  in  H^(Ri  U  Kj)  C  H^fR,)  x  H^{Kj),  extended  by  zero  outside.  Then  the 
boundary  terms  J{u,v)  and  J{v,  u)  vanish,  because  [m]  =  [u]  =  0  on  jy,  and  the  weak 
formulation  (3.16)  reduces  to 


Vv  +  cuv)dx=  /  fvdx 

JKiUKj 


(3.30) 


On  the  other  hand,  multiplying  (3.29)  by  v,  integrating  on  Kj  and  Kj  and  using 
Green's  formula,  we  have: 


f  CVu-'Vv  +  cuv)dx-  [  (n  ■  V u.)iV ds  =  /  fvdx, 

JKi  J-iii 

f  (Vu-'Vv  +  cuv)dx—  I  {n-Vu)jVds  —  /  fvdx, 
JKj  Jlii 


so  that 


f  (Vu-Vv  +  cuv)dx-  f  [n-VM]i;ds=/  fvdx. 

JKiUKj  ha 


(3.31) 


Comparing  (3.30)  and  (3.31),  one  observes  that: 

/  [n-VM]t;ds  =  0,  Vu  G  H^(Ki  U  Xy). 

J'lii 

Then,  [n  ■  Vm]  =  0  for  all  element  edges  jy,  which  implies  Vm  G  H{div,Q).  This 
allows  us  to  conclude  that  u  satisfies  Poisson  Equation  globally  on  Q,  i.e. 


-Am  +  cu=f,  a.e.  in  Q. 


(3.32) 


To  recover  the  Dirichlet  boundary  conditions,  we  now  consider  a  function  v  G 
H^(a)  n  H2(Q),  so  tiiat  integrating  (3.32)  provides: 


f(Vu-Vv  +  cuv)dx=  fvdx, 

Jci 


whereas  (3.16)  yields: 


f  (Vu-'Vv  +  cuv)dx-  f  in-Vv)uds=  [  fvdx-  (n-Vv)uods. 

Ja  Jfd  J^D 


Substracting  both  equations,  we  obtain: 


f  (n- Vt;)(M-Mo)ds  =  0,  \/v  e 

JTn 


and  conclude  that  m  =  mq  on  Tp. 

In  the  same  way,  choosing  v  G  C  such  that  u  =  0  on  Tp,  we  get: 

f  {n-Vu- g)vds  =  0, 

Jtn 

sothatn- Vm  =  gonrN.  G 

Remark  4  When  c  is  zero,  C^(0)  can  he  replaced  in  Theorem  3.1  by  D  H^i^h)  since 
Vu€H{div,Q). 
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3.4.  Properties  of  the  Bilinear  Forms 

3.4.1.  Mesh-dependent  norms 

We  now  introduce  norms  associated  with  the  bilinear  forms: 

1.  Energy  Norm: 

Ml(P„  =  'B{v,v)=  X  ll^ll6,K=  X  (|I^^IIo,k  +  ‘^II^’IIo,k)  (3-33) 

KeSPfc  Ke55,  ^ 

2.  Norm  proposed  by  Siili  et  al.  in  [22,15]: 

\\v\\\  =  B{v,  v)  +  fiv,  v)  =  ||z;||2  ^  +  f  (T  [vf  ds  (3.34) 

•'Il'nfUrD 

3.  Norm  proposed  by  Baumann  et  al.in  [7,17,18]  and  by  Baker  and  Karakashian 
in  [6]: 

M\k=Mk+L  ,  Ut>-Vvf<k  (3.35) 

•'IfwfUlD  ^ 

We  note  that  the  energy  norm  becomes  a  seminorm  when  c  is  zero. 

3.4.2.  Continuity  of  the  bilinear  forms 

We  shall  show  now  that  the  bilinear  forms  ®i(-,  •)  and  f^(-,  •)  are  continuous  on 
H\%)  with  respect  to  the  norm  IIHIlajj^  defined  in  (3.35).  Unfortunately,  we  are 
tonable  to  show  continuity  with  respect  to  the  other  two  norms  (3.33)  and  (3.34). 

Theorem  3.2  (GEM  and  DGM)  Let  ‘S±{-,  •)  he  the  bilinear  form  defined  either  in  (3.15) 
or  in  (3.23)  .  Then, 

mu,v)\  <  \\\u\\\^^  |||z;|||25, ,  Vm.Vi;  €  H^tP,,).  (3.36) 

Proof: 

First  note  that: 

\'B±{u,v)\  =  \B(u,v)-J{u,v)±J(v,  u)\ 

<  \B(u,v)\  +  \J{u,v)\  +  \Jiv,u)\ 

It  is  clear  that 

\B{u,v)\<  X  [  IVm- Vu  +  cwyldr  <  ||w||e,2^J|u||e,2}, 

KeS), 
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The  first  boundary  term  gives: 

\J{u,v)\<  f  |(n- VM)[z;]|ds 


<  ,fj  0-1  (n-V«)^ds  ./7  (X  [t;]^ds. 

“  V  -/n«<urD  V  JrMoro 

Likewise, 

i(n-Vt;)[M]|ds 


<  ,  /7  <T  [wf  ds  ,  f[  <T-'^  (n-Vt;)^  ds. 

V  Jr-Murp  V  -/i/njurD 

In  consequence,  we  have,  using  the  discrete  Schwarz  inequality  (A.l): 
\(B±{U,V)\  <  \\u\\e,(B,\Me,% 

+  .fj  0-1  (n-VM)^ds  ,/7  <7  [o]^ds 

V  Jr^^Tp  V  •'ifnfUrD 

+  *[7  o-  [M]^ds  ,/7  0-1  (n-Vo)^ds. 

V  ^ri„,urD  V  ‘/TlnfUrD 


X 


<  IWU  lli-lls , 


which  completes  the  proof. 


□ 


Theorem  3.3  (SIPG  and  NIPG  Methods)  Let  •)  be  the  bilinear  form  defined  ei¬ 
ther  in  (3.18)  or  in  (3.26).  Then, 

l®l(".o)l  <  Cll-lllai  iKlk  .  Vu.Vi.  e  H\%).  (3.37) 

where  C  is  a  constant,  C  <2. 

Proof: 

As  before  we  have: 

v)\  =  \B{u,  V)  -  }{u,  V)  ±  }{v,  u)  +  r{u,  0)1 
<  |B(m,o)|  +  \Ku,v)\  +  1/(0,  «)|  +  |/<^(u,o)| 

<iiiM|ikiMi2>,  +  in».^)i- 
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And 


ir(M,  1^)1  <  f  W  [u]  [uJI  ds  <  .fj  a[ufds  .fj  ds. 

•'riniuro  V  •'riniuro  V  •'i/nfUrD 

Therefore,  making  use  again  of  the  discrete  Schwarz  inequality  (A.l),  we  obtain: 
\‘^iu,v)\  <  IuIt^  +  \l L  ^  (T[ufds  J  f  a[vfds. 


^  2  lllwllly^  > 

and  we  see  that  C  is  at  most  equal  to  2.  □ 

3.4.3.  Coercivity  of  the  bilinear  forms  in  the  discrete  spaces 

Here  we  wish  to  show  that  the  bilinear  forms  ®t(-,  •)  and  •)  are  coercive  in 
H^((Ph)  with  respect  to  the  norm  |||-  |||5j_  in  order  to  be  able  to  apply  classical  theorems 
for  existence  and  uniqueness  of  solutions  of  the  discontinuous  methods.  Unfortu¬ 
nately,  to  date,  we  are  able  to  prove  coercivity  only  in  the  discrete  discontinuous 
spaces  and  then,  only  for  the  SIPG  and  NIPG  formulations  . 


Theorem  3.4  (NIPG  Method)  Let  a  =  np^/h,  k  being  a  positive  number.  Then,  for  all 
K>  0,  there  exists  a  positive  constant,  a>0,  such  that: 

%{z,  z)>oc  |||2||||_ ,  Vz  €  (3.38) 

Here  a  is  independent  ofh  and  p. 


Proof:  Let  a  be  an  arbitrary  real  number  and  choose  a  z  e  Then 
S;(z,z)-a||z||j 

=  (1  -  a)  B(z,z)  +  (1  -  a)  J'’{z,z)  —  af  —  (n  ■  V2)^ds 

./r-„,urD  cr 

Since  (n  •  Vz)  is  the  average  of  the  flux  at  the  interface  of  two  elements  K,-  and 
Kj,  the  corresponding  integral  can  be  split  into  two  integrals  with  integrands  (n  • 
Vz)ija  and  (n  •  Vz)j/a,  each  one  associated  with  the  elements  K,  or  Kj  respectively. 
Therefore,  let  7  ^  ^int  U  Td  and  consider  the  integral  associated  with  the  element 
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K.  Using  the  trace  inequality  (A.3)  and  the  inverse  property  (A.7),  we  have 

f  -(n- Vzf  ds  <  ^||Vz||g 
J-y  or  <T 

<§(i||Vz||?,K+||Vz||„,,:||V2z||o,K) 

so  that,  selecting  a  to  be  equal  to  K'p^/h^,  we  obtain: 

-/i(iiVzfds>-2||Vz||J,K. 

Jt(t  k  ’ 

Note  that,  when  the  mesh  size  and  /ik;  and  the  polynomial  degrees  Pk,.  and 
are  different  from  each  other  in  the  two  elements  K,  and  Kj  sharing  the  edge  7,y,  we 
actually  choose  cr  as 

max(p|.,p|p 
^  ~^min(/ij;., /ijc.)’ 


so  that: 


2 

0,K, 


K  max(p|.,  hKi 


2 

0,Ki 


<  §l|Vz||? 


, Ki¬ 


lt  then  follows  that: 

^^{z,  z)-a  \\\z\\\l^  >(l-a-  aC/n)  B{z,  z)  +  (1  -  a)  T  (z,  z). 
Therefore,  we  certainly  can  pick  a  value  of  a  such  that 


for  which  the  bilinear  form  •)  is  coercive  in  for  all  «  >  0.  □ 
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Theorem  3.5  (SIPG  Method)  Let  a  =  Kp'^/h,  n  being  a  positive  number.  Then,  for 
K  >  kq,  there  exists  a  positive  constant  a  independent  ofh  and  p,  a  >  0,  such  that: 


Proof:  Let  a  be  an  arbitrary  real  number  and  choose  z  G  Then 
z)-a  \\z\C  =(!-«)  B(z,  z)  +  (1  -  a)  r(z,  z) 


(3.39) 


-if  (n  •  Vz)  [z]ds  —  af  —  (n  •  Vz)^ 
^n„,urn  (T 


ds 


«fUrjc>  JTintUVo 

There  exists  a  positive  number  e  such  that  for  every  edge  7  €  lint  U  Fd: 


2 (n  -  Vz)  [z] ds  <2  ^ J  {n^Vz)^ds^J J  (7[z]^ds 

<e  /  —  (n  •  Vz)^ds  +  -  f  a[z^ds 
J-y  a  e  J'y 

which  yields,  using  the  result  in  the  previous  proof: 

(z,  z)  -  a  |||z|||y^  >  ^1  —  a  -  (a  +  ^)~)  ® 


In  order  to  prove  coercivity,  we  want  to  find  o  >  0  such  that  both  factors  in  the 
inequality  are  positive,  in  other  words: 


—  a  —  (a  +  ^)~)  >  0  and  —  a  —  >0. 


The  second  inequality  requires  that: 


0<a<l-- 
“  e 

which  means  that 

e  >  1. 

On  the  other  hand  the  first  inequality  requires  that: 

„  ^  \  —  eC! K  _  1  —  C/ K  .  K  —  C 

l  +  C/«  ~l  +  C/n~K  +  C 

This  completes  tf\e  proof  by  taking  «  sufficiently  large,  namely  k>  kq  (where  for 
instance  «o  >  C.)  □ 
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Remark  5  We  note  that  •)  (for  NIPG  Method)  is  coercive  in  H\%)  with  respect  to 
the  norm  H-Ha},.  Indeed,  for  all  v  € 

>Sf.{v,  v)  =  Biv,  v)  -  J(v,  v)  +  Kv,  v)  +  r{v,  v)  =  \\v\\\ .  (3.40) 

It  is  also  straightforward  to  show  that  ^+{- ,  •)  (for  DGM)  is  coercive  in  H^ifPh)  with  respect 
to  the  energy  norm  IHIe.a),- 

<B+{v,  v)  =  B(v,  v)  -  J{v,  v)  +  J(v,  v)  =  B(v,  v)  =  \\vfe,% ■  (3-41) 

These  results  will  be  crucial  in  deriving  a  priori  error  estimates  in  the  next  section. 


4.  A  Priori  Error  Estimates 


4.1.  SIPG  and  NIPG  Methods 

Theorem  4.1  Let  u  e  H^(Q)  D  IP(%),  s>2,bea  solution  of  (3.18)  (SIPG)  or  (3.26) 
(NIPG)  and  uh  be  the  discrete  discontinuous  solution  of 


'BliUh,v)=m,  (4.1) 

Then,  choosing  a  -  np^fh,  (k  >  O/or  NIPG  and  k>ko  for  SIPG),  the  numerical  error 
e  =  u-Uh  satisfies: 


where  p  =  min(p  +  l,s)  and  p  >  1. 


(4.2) 


4.1.1.  Proof  of  Theorem  4.1  for  SIPG  and  NIPG 

First,  by  definition  of  the  norms,  we  note  that  ||e||e,25,  <  |||e|||^.  In  other  words,  it 
suffices  here  to  estimate  the  error  with  respect  to  the  norm  Ill-IHy^.  The  proof  is 
inspired  by  [5,6,16]  where  the  authors  have  derived  the  rate  of  convergence  in  h 
only  for  the  SIPG  method  of  the  (3.22)  form.  Here  we  extend  their  results  to  the 
NIPG  formulation  as  well  and  also  show  for  both  methods  the  rate  of  convergence 
in  p. 

Proof:  Let  Zp  be  an  interpolant  of  u  in  We  shall  use  the  notation  rj  =  u  —  Zp 
and  ^  =  Uh-Zp  so  that  e  =  u-Uh^r]-^.  Applying  the  triangle  inequality,  we 
have: 

Mk  =  ill«  -  ^h\h  =  111^  -  ^\h  ^  • 
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From  the  coercivity  of  the  bilinear  form  ®±(*,  •)/  since  ^  G  we  have 

IKlis  < 

and  from  the  "orthogonality"  property  —  Uh,v)  =  0,'iv  e  we  get 
Using  tire  continuity  of  ®±(-,  •)/  we  know  that 
which  implies 


Kllk<C|||^il2.,. 


Finally,  we  have 

IWl3i<Wk  +  lll«llk<c||l-)llj|,. 

We  recall  here  that  C  is  a  generic  constant  independent  of  h  and  p  which  takes 
different  values  at  different  places. 

We  now  choose  the  interpolant  Zp  as  defined  in  Lemma  A.7.  Then: 


lll’)llii=  I  X(|V,p+c,^)dA:+/ 

Ke'Pk''^ 


ifuro 


—  (n-V7j)^ds+  f  (7[77]^ds 

(4.3) 


The  integrals  in  the  leading  term  are  estimated  as,  using  (A.8): 
^|V,pdx<c(^^  ||u|p,,.  s>l, 

l|i<lP,i:.  s>0. 

SO  that 

J^(\Vr]\'^  +  cr]^^dx<C^^\\u\\lK,  s>l. 

^  P  K 

Let  ')ij  denote  an  interior  edge  shared  by  the  elements  Ki  and  Kj.  Then,  using  the 
inequality  {a  +  Vf  <  2a^  +  we  observe  that 

[  -{rx-  Vrjf  ds<l  [  -  (n  ■  ds  +  i  ^  (n  •  (Vr/)  )  "  ds. 

^7,;  2  Jjij  a  I  Jyij  a  \  ‘J 
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In  other  words,  in  splitting  the  second  integrals  on  the  right  hand  side  of  (4.3)  as 
above,  we  actually  associate  with  each  7  €  'Eh^int  element  K,  such  that 


J-v  a  cr 


c  ( 1  Ifr^  .  K-' 


C  ( II  1,2 

<  -  -^=2  +  “fe  “  ls,K 

\Pk  Vk  J 


C'  ^ 


s>2. 


Again,  for  an  interior  edge  7,;  shared  by  K,-  and  Kj,  using  {a  -  hf  <  2a^  +  2}p-,  we 
have: 

f  a[rtfds^f  (r{rii-'njfds<2  f  (7{'nifds  +  2f  a  {r}jf  ds 

J^ij  hij  J'iii 

This  means  that  the  edge  integrals  making  the  third  term  of  (4.3)  are  bounded  by: 

f  (r(r]fds  <  C<t^^\\u\\Ik  <  C«-fc||M||2^ 

•'7  Pk  Pk 

In  combining  the  above  results,  we  thus  obtain 


r  ^2p-2  ^2/i-2  ^2ix-2 

i  C  I  -feTillBlUK 

Keff,  Pjc 

s  c^ll"ll» 

which  is  the  expected  a  priori  error  estimate. 


-2  'k  1/2 
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4.1.2.  Alternative  Proof  of  Theorem  4.1  for  NIPG 

Alternatively,  we  present  a  second  proof  of  Theorem  4.1  for  tiie  NIPG  method  only 
as  it  is  based  on  the  nonsymmetry  of  the  formulation.  The  proof  is  inspired  by 
the  one  found  in  [22].  However,  our  rate  of  convergence  with  respect  to  p  was 
improved  from  (s  —  2)  to  (s  -  3/2)  using  the  interpolation  estimates  of  Lemma  A.7. 
Later,  the  same  authors  proposed  in  [15]  a  comparable  version  of  the  proof  with 
(s  —  3/2)  as  the  rate  of  convergence. 

Proof:  Once  again,  Zp  is  the  interpolant  of  if  in  as  defined  in  Lemma  A.7.  and 
we  denote  ri  =  u  —  Zp  and  ^  —  Zp  as  before.  Then, 

<  Ikikft  =  II w  -  W/illPfc  =  Ih  -  ^1125,  <  MWh  +  11^1125, • 

Moreover,  from  the  definition  of  •)  and  the  norm  IHIe.a),  (see  (3.40))  and  the 
"orthogonality"  relation,  we  have: 


The  goal  is  now  to  bound  ®+(r7, 0  ii'  terms  of  .  Recall  that: 

mv,  0  =  Bin,  0  +  nn,0-  Jin,  0  +  /(C,  n) 

<  \Bin,o\ + + \jiri,o\ + m,v)\ 

The  first  term  on  the  right  hand  side  of  the  equation  above  gives: 

|b(7?,OI<  X  f  \^n-^^  +  cn^\dx<\\n\\e,%\me,^,<\\nU^^^ 

The  term  iJ'^in,  0|  is  bounded  by: 

ir(77,ai<  /  kMmids 

JTiniUTD 

<  xfl  <x[nfds.fj  <r[^]^ds 

V  Jr^urD  V  JrintUTD 

<  ll^lkjl^lk*. 

whereas  we  have  for  the  third  term: 

\m,0\<  f  Kn  V,,)Kl|d5 
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Likewise,  /(^,  rf)  is  bounded  by: 

|/(^,^)|  <  MkJ  L  ^  <T-i(n-VO^ds 

V  ''WnfUFo 

Using  again  the  trace  inequality  (A.3)  and  the  inverse  property  (A.7),  it  is  shown 
that: 


In  other  words,  using  <t  =  Kpj^/ht: 

|/(^,T/)|  <  C||r/||!pJ|C||2), 

Combining  the  above  results,  we  have: 

IlClk  ^  ^  ^ll^lk  +  ^  j  ^ 

so  that: 

Iklle.y*  <  IkIliPfc  <  ^  Ibik  +C|||»7|is5,  -  • 

We  conclude  the  proof  by  employing  the  estimate  on  |1|//|||2,^  shown  in  the  previous 
proof. 

4.2.  DG  Method 

We  recall  that  the  DG  formulation  proposed  in  [7,18]  is  deduced  from  the  NIPG 
method  by  simply  setting  the  penalty  parameter  o  to  zero.  However,  unlike  NIPG, 
continuity  and  coerdvity  of  the  bilinear  form  ®+(-,  •)  cannot  be  proved  simultane¬ 
ously  using  the  same  norm.  At  best  it  is  shown  that: 

and  ti\at: 

«+(».*>)  Ms,.  v«,!,6hW 

The  main  issue  in  finding  a  priori  error  estimates  for  the  error  e  =  m  -  «/,  in  the  nu¬ 
merical  approximation  of  the  DG  problem  consists  in  deriving  an  upper  bound 
on: 
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with  respect  to  the  norm  whenc  =  0.  This  integral  does  indeed  appear  when 

bounding  the  term  J(t],  0/  i-e. 

\I(V.()\<  f  |(n-V,,>({J|ds<.//;  {«  VrifisJf  Kfds. 


We  present  below  two  approaches,  by  treating  separately  the  case  when  c  is  zero 
and  the  case  when  c  is  nonzero. 

4.2.1.  A  priori  error  estimate  when  c  is  nonzero 

We  find  it  instructive  to  analyze  the  special  case  in  which  c  is  strictly  greater  than 
zero.  In  this  case,  we  still  can  use  the  methodology  presented  earlier  for  the  NIPG 
method.  However,  we  shall  see  that  the  rate  of  convergence  with  respect  to  the 
mesh  size  becomes  suboptimal  as  stated  in  the  following  theorem. 


Theorem  4.2  Let  u  e  H^(Q)  n  s>2,bea  solution  of  (3.23)  with  c  >  0  and  he 

the  discrete  discontinuous  solution  of  (3.24).  Then,  the  numerical  error  e  =  u  —  Uf,  satisfies: 


<  C 


h^^ 


-2 


'\e,!Ph  ^  ''jjS-3/2' 


(4.4) 


where  p  =  min(p  + 1 ,  s)  and  p  >  1. 


Proof:  Using  the  same  procedure  and  notation  as  before,  we  have: 

Iklle.Sj,  =  I|m  -  Uh\\e,!P,  =  11^  “  i\\e,%  <  h\\e,V,  + 

Moreover,  from  the  definition  of  ®+(-,  •)  (see  (3.41))  and  the  "orthogonality"  rela¬ 
tion,  we  further  show  that: 

=  ®+(^,a 

<\B(v,0\  +  \J(il,0\  +  \J(Ori)\ 

We  now  consider  each  term  one  at  a  time.  The  first  term  B(r],  ^  is  straightforwardly 
bounded  by: 

|B(7y,a|  <  \\v\\eMeM<C-^\W\\M\\e,n 


(4.5) 
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We  expect  that  the  third  term  J{^,t})  can  be  treated  as  before  and  should  not  pose 
any  problems.  Indeed,  applying  the  Cauchy-Schwartz  inequality,  we  have: 


When  7  C  T^t  U  To  and  C  ^  we  have  already  shown  tiiat: 

l^{n-VO^ds<C^m\\lK- 

Next,  we  obtain  from  the  approximation  property  (A.9) 

f  r]^ds  =  \\r]\\l^  <  C-^||w||^,K- 

Pk 

Therefore  the  term  J(^,  rj)  is  bounded  by: 


Finally  we  need  to  consider  the  term  which  is  held  responsible  for  deterio¬ 

rating  the  convergence  rate  of  the  solution.  By  the  Cauchy-Schwarz  inequality,  we 
have: 

Once  again,  the  approximation  property  gives 

/(n.V;,fds<C^||«||^.K 
Pk 

while  from  the  trace  inequality  (A.3),  we  have: 


l«ll5.,<c(i||«iiS,K+iieiio,,iiv«iio,K) 

<c(i||{|l§,,+i|ieilS,K+/.Kilv{|i5,,) 
<c(^ii«iiS,k+''kIIv«ii5,,) 
<c(i|l?ll5,,+ /Idling, k) 

<  i-mi. 


(4.6) 
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It  is  important  to  point  out  here  that  the  norm  |l^i|o,7  is  bounded  as  long  as  c  >  0. 
Then  we  have: 


h' 


In  conclusion, 


ll^lle.25-  ^  ^  +  31375  +  3372 


ua-2 


which  completes  the  proof. 


□ 


Remark  6  Note  that  C  is  inversely  proportional  to  c.  Therefore  the  error  is  expected  to 
grow  as  c  gets  smaller. 


4.2.2.  Discussion  of  the  case  in  which  c  is  zero 

The  operator,  when  c  is  zero,  reduces  to  the  pure  Laplacian.  In  this  case,  the  energy 
norm  \\-\\e,%  becomes  the  seminorm  ||V-||o,!P*-  Following  the  same  procedure  as 
before,  we  would  have: 

l|V^||g,2>,  =  ®f(e,0  =  ®+(»7,0  (4.7) 

where  t]  =  u  —  Zp,^  =  Uh  —  Zp  and  Zp  defines  an  arbitrary  interpolant  of  u  on 
However,  from  (4.6),  we  can  see  right  now  that  the  term  ®+(7/,  0  would  then  be 
bounded  by  l|^||o,a5,.  In  turn,  it  is  impossible  to  bound  ||^||o,2’^  with  respect  to 
||V^||o,25,.  Therefore,  the  previous  methodology  to  obtain  error  estimates  cannot 
be  applied  in  the  present  case. 

Suppose  that  we  introduce  an  elementwise  constant  function  ^  to  be  defined  later. 
Then,  we  can  rewrite  (4.7)  as: 

II  V^loV.  =  0  =  'S+(rj,  e  - 1  +  0  =  ®+(77,  e  -  0  +  'S+iri,  0-  (4-8) 

Suppose  now  we  can  construct  a  new  interpolant  such  that: 

'S+(v,O^0.  (4.9) 


Then  we  would  have 


(4.10) 
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We  have  seen  that  the  terms  5(7?, 6  and  J{^,ri)  are  easily  bounded  in  terms  of 
II  V^||o,25,-  The  other  term  reads: 

}iV,^-0=  f  <*fVr/)[^-|]ds. 

According  to  Lemma  A.5,  this  integral  can  be  bounded  with  respect  to  ||V^||o,35, 
under  the  condition  that  ^  is  chosen  as  the  average  of  ^  on  each  element. 


This  approach  has  been  followed  in  principle  by  Riviere,  Wheeler  and  Girault 
in  [20,19]  where  they  construct  special  interpolants  ttm  which  satisfied  (4.9)  and 

||M-7rM||o,K<C-;^||Ml|s,K, 

Vk 


||V(m  -  7rw)||o,K  <  C-fi2 
Vk 


||V^(m  -  7rM)||o,fC  <  C-fi2  ll“lls,K! 

Vk 

where  /z  =  min(pK  +  1,  s),  s  >  2,  pK  >  2.  Using  these  interpolants,  they  were  able 
to  derive  an  a  priori  error  estimate  of  the  form: 


||Vel|o,2i<C^||M||5. 


(4.11) 


Although  the  rate  of  convergence  is  optimal  in  h,  we  show  next  that  the  rate  of 
convergence  in  p  is  in  reality  better  than  (s  -  4).  We  improve  this  result  by  con¬ 
structing  better  approximation  properties  for  the  new  interpolant  and  by  refining 
the  analysis. 


4.2.3.  New  Interpolants 

Lemma  4.1  Let  Kbea  triangular  element  of  the  partition  tPh  and  u  a  function  in  H%K), 
s  >  2.  There  exists  a  positive  constant  C  depending  on  s  and  g  but  independent  of  u,  Pk, 
and  hji,  and  a  polynomial  ttm  G  Pp^^iK),  px  >  2,  such  that 


and 


f  n  •  V(m  —  ttm)  ds  =  0,  V7  C 

•'7 


h^ 

\\u  -  7rM||o,K  <  c  s_3/2ll“lls,K) 

Vk 

||V(m  -  7rw)||o,K  <  C-^§3y2 
Vk 

||V^(m  -  7rM)||o,K  <  C-^\\u\\s,K, 
Vk 


dK, 


(4.12) 


(4.13) 
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Figure  2.  Reference  element  K  and  mapping  Fk  from  K  to  the  ele¬ 
ment  K  in  the  physical  domain. 


where  n  =  min(pjc  + 1>  s). 


We  present  the  proof  of  this  theorem  for  triangular  elements  only.  The  proof  is 
similar  for  quadrilaterals. 

Proof:  Let  the  triangle  fC  €  Q  be  the  image  of  the  master  element  by  the  affine 
mapping  Fk  as  shown  in  Figure  2.  The  mapping  is  often  rewritten  as: 

FK(x)  =  Bx  +  b  (4.14) 

where  B  represents  a  two-by-two  matrix  whose  components  are  independent  of  x 
and  b  is  a  two-dimensional  vector.  Here,  7  will  refer  to  the  edge  between  node 
N2  and  N3,  unless  stated  otherwise,  and  7  on  K  will  denote  its  image  by  We 
associate  with  7  and  7  the  unit  normal  vector  n  and  n,  respectively. 


Given  rj  6  H^(fC),  namely  rj  =  u  —  Zp,  where  Zp  is  the  interpolant  of  u  as  defined 
in  Lemma  A.7,  the  objective  here  is  to  construct  a  polynomial  function  q  in 
such  that: 


/  n  •  V?/ds  =  /  n  •  Vqds. 
t/7  J'y 


(4.15) 


Indeed  we  would  have: 


/  n  •  V{ti  -q)ds=  /  n  •  V(u  —  Zp  —  q)ds  =  /  n  •  V(m  -  {zp  +  q))  ds  =  0, 

77  J'y 


and  the  new  interpolant  could  be  derived  as  ttm  =  Zp  -f  q. 


Following  [19],  and  assuming  ^  2,  we  introduce  the  polynomial  function  q-y 
associated  with  the  edge  7  on  K: 


qy  =  C^(l  -X-  y){x  -I- y), 


Vx^{x,y)ek. 


(4.16) 


Figure  3.  Polynomial  function  on  the  reference  element  K. 

where  C-y  is  a  constant  to  be  defined.  We  observe  in  Fig.  3  diat  such  a  polynomial 
function  satisfies: 


/  fi  ■  ds  =  0 

>'731 

with  7,y  defining  an  edge  on  K  joining  the  nodes  N,-  and  Nj. 

The  constant  is  found  so  that  (4.15)  is  satisfied  in  the  physical  space.  Further¬ 
more,  we  obtain  an  upper  botmd  for  Cy  (see  Riviere  [19])  as: 

|C,|  <  C||Bf  IIB-'f  IIV^IIo,, 

<C<r^0)  /if  l|V.)||o,, 

<C/if  IIVbIIo,, 

<  C  |||V77||o,k  +  ^kI|V??||o,kI|V^^||o,k} 
where  we  make  use  of  the  Trace  Inequality  (A.3).  We  also  observe  that: 

\\q4o,K<C\detBnqy\\o,^ 

<ChK\Cy\\K\ 

<  Ch}^\C'y\> 
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Likewise,  we  have: 
liV<?,||o,K  <  C|C,| 


So  far,  we  have  carried  out  the  analysis  for  the  edge  7  between  node  N2  and  N3. 
We  point  out  that  the  same  results  are  obtained  for  the  other  two  edges.  We  then 
associate  with  each  edge  712,  723/  731  /  a  polynomial  ^12,  ^31  respectively  such 

that 

^i2  =  Ci2y(l-y) 

=  C23  (1  —  x  —  y){x  +  y) 

^31  =  Qi  ^(1  - 

Adding  these  polynomial  functions  together,  we  construct  on  the  element  K  a  new 
function  q  €  PiiK) 

<7(x)  =  qni^)  +  '?23(x)  +  <?3i(x),  Vx  €  JC, 
which  satisfies 


/  n-'Vqds=  /  n  ■  V((^i2 +  (/23 +  ?3i)ds  =  /  n-Vq^ds^  /  nVrjds, 

‘'712  *'712  •'712  ''712 

/  nVqds=  n  •  V(^12 +  <?23 +  ^3i)ds  =  /  n-'Vq23ds=  n-Vqds, 

>'723  •'723  •'723  •'723 

/  n-Vi?ds=  /  n-V(«ji2  +  <723  +  ^3i)ds=  /  n-V^3ids=/  n-Vr/ds. 

‘'7.11  ‘'731  ‘'731  *'731 
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Now,  by  the  triangle  inequality, 

||m  -  TTu\\o,K  <  \W  -  Zp||o,K  +  Ikllo.K 

<  Iloilo, K  +  lkl2||o,K  +  ll<?23||o,JiC  +  Iksillo.K 


< 


||»?||o,K  +  C/zk  |||Vjy|lo,K  +  1jj<:||V?7||o,k||  V^r?||o,K  j 


1/2 


<C 


ht 


K 


'hl^ 


Pk'^^^KPk  ^'^^^Pk^Pk^ 


}  l|M|kK 


h'^ 

<C-^ 


H-1  ) 

hfi  I  ll“IU 


wHs.ic- 


Pk 

In  the  same  manner,  we  find: 

||V(m  -  7rM)||o,K  <  II V(m  —  Zp)\\o,K  +  ||V(7||o,K 


<  II V7y||o,K  +  c  {|| V77||g,K  +  /»kII V;?||o,k|| V^T/lkx} 


1/2 


^  K 

[Pk 


+ 


pr'\ 


«lkK 


<  c 


-1 


P 


iS-3/2 


|W|Is,K5 


and 


II  V^(M  -  7rM)||o,K  <  II  V^(M  -  Zp)||o,K  +  II  V^^||o,K 


<  II  V^»7||o,K  +  Chf}  I II  Vt^IIo^ji;;  +  /»kII  V?7||o,kII  V^»7||o,k} 


1/2 


^  f  K~^  ,  1  K~^  1 

<  C^Ms,k- 


m\s,K 


We  observe  that  the  first  two  estimates  are  governed  by  the  rate  of  convergence  of 
WqWo^K  and  II  V^7||o,K  respectively,  while  the  last  estimate  is  governed  by  the  rate  of 
convergence  of  ||V^77||o,k*  D 
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4.2.4.  A  priori  error  estimate  when  c  is  zero 

Theorem  4.3  Let  u  €  D  s>2bea  solution  of  (3.23)  and  Uh  be  the  discrete 

discontinuous  solution  of  (3.24)  with  c  =  0  and  p  >2.  Then,  the  numerical  error  e  = 
u  —  Uf,  satisfies: 

(4.17) 

where  p  =  min(p  +  1,  s). 

Proof:  Let  ttw  be  the  interpolant  of  u  in  defined  on  each  element  fC  of  2^  as 
in  Lemma  4.1.  We  also  introduce  ri  =  u  —  iru  and  ^  =  M/i  —  ttm.  Using  the  triangle 
inequality,  we  have: 

II Velio, =  II V(77  -  OII0.25,  <  II V77||o,25,  +  II VCHo.iP, 

and  from  (4.7)  and  (4.8),  we  recall  that: 

II  V^||g,2>,  =  ®+(^,  0  =  ®+(»7,  ^  -  0  +  'S+M- 

Here  |  is  chosen  as  the  average  of  ^  over  each  K,  i.e. 


We  note  here  that  the  authors  in  [20,19]  chose  |  as  the  average  of  ^  over  each  edge 
and  their  proof  is  thus  slightly  different  from  ours. 

This  particular  choice  of  the  interpolant  ttm  and  piecewise  constant  function  ^  does 
indeed  yield; 

'B+ip,  0  =  B(rj,  0  -  J(r),  0  +  }(l  rj)  = -J(v,  0 

=  -/  (n.V7?)[e]ds 

JTmUTd 

=  -[^]  f  (n-V7/)ds 

•'Ij'ntUro 

=  0 

since  the  last  integral  is  zero  according  to  the  property  (4.12)  of  the  interpolant  ttm. 
Therefore 

l|V{|l5,s,.  =  «+(>!.e-f)  (418) 

We  now  show  how  B+(r],  ^  —  f)  can  be  bounded  with  respect  to  ||  V^||o,25,.  We  nat¬ 
urally  have  from  (4.10) 

i®+(r/,e  - 1)1  <  iB(^,oi + \jiv,^- 1)1 + \mv)\ 
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The  first  term  gives,  using  the  approximation  properties  of  Lemma  4.1  and  the 
discrete  Schwarz  inequality: 


X  f  |V77-V^|da:<  ^  ||V7?||o,kI|V^||o,k 


Ke®  Pk 


The  third  term  is  treated  as  usual.  We  have 

I7(e,r?)|  <  /  Kn- Ve)[/?]|ds<XlKn- VOIIo,7llMllo,7 

<CX  X  lln-V^||o,7l|r?||o.7 

K€<p^')€dK\rtt 

<cX  X  llveilo,7Nlo,7 

K€2’fc7€9K\rN 

From  the  trace  inequality  (A.3)  and  the  inverse  property  (A.7),  we  show  that: 

r  1  1 

<  c|i||V{|g,,+  ||V{||„,,C„^||V«||„,,}  ' 

<c^mhK 

and,  from  the  approximation  properties  of  Lemma  4.1: 

r  1  1 

lkllo,7  <  c|^lhllo,f:  +  II»7||o,k||V»71|o,j:c| 


<c/^L  +  ^ 

-  ^  1  „2s-3  „2s-3 

(Pk  Pk 
^,^*-1/2 


Pk  '  Pk 
-1  1/2 


Plkfc 
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In  conclusion,  we  find  that: 

l/({.>!)l<c  2  ;^l|V«llo,K  ^11 

K€2), 

^  C  2  ^IIV«llo,dl»lkic 

KefSipK 


w||s,K 


In  the  same  manner  as  before,  we  obtain  for  the  term  ^  —  ^)\ 


|/(r/,C-|)|<C  X  S  l|Vr/||o.7ll^-|||o,7 

K-e2),  -i€dK\Vt^ 

In  this  case,  we  have  using  also  the  approximation  properties  of  the  interpolant: 

1/2 


l|V.,||c,,  <  C  {  ^l|Vr,|g,K  +  II  V,||o,k||  V"^||„,k} 

f  1  „2  II  Il2 

-  ^s^-3/2  ^^-2  il"lkK 


1/2 


<c 


_ , 

pJ-3  + 


}'wi. 


;,M-3/2 

Vk 


However,  for  the  other  term,  we  have,  using  Lemma  A.5: 

Ilf  -  f  llo,,  5  C  |L||f  -  f  ||5,k  +  Ilf  -  f  ||o,i!||VK  -  f)||„,it| 

31/3 


1/2 


<  l^llf  -  f  llo,f:  +  Ilf  "  f  llo,K!|Vf  llo,K  j- 

f  1 

<  C  |i).i:l|Ve||J,K  +  f.d|Vf||o,d|Vf  I|o,k| 

<  af  liVf  ||o,K 
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It  follows  that: 


;,M-3/2 

I7(>),«-?)I<C  S 

K€!PhPk 

<  C  I  ^ll«lkitl|V«l|o,i; 
5  C^I|a|W|V«||„,:,l 


Combining  the  previous  results,  we  finally  get: 


pS-7/4  +  ^-5/2  ) 


l|w|ls<C 


pS—5/2 


ll«lls 


and  this  completes  the  proof  since  ||  V?7|lo,!Pj  converges  with  a  greater  rate  of  con¬ 
vergence  than  ||V^||o,Pj,-  Cl 


4.2.5.  Alternative  estimate  when  c  is  nonzero 

We  now  use  the  previous  results  to  review  the  error  estimate  when  c  is  nonzero. 
The  new  estimate  is  given  in  the  following  theorem: 


Theorem  4.4  Let  u  G  n  s>2,bea  solution  of  (3.23)  with  oO  and  he 

the  discrete  discontinuous  solution  of  (3.24).  Then,  the  numerical  error  e  =  u-Uh  satisfies: 


where  //  =  min(p  -1- 1,  s)  and  p>2. 


(4.19) 


Proof:  In  this  case,  we  have: 

=  B{rj,0-J(V,0  +  mr}) 

=  Biv,  0-J(r),C-0-  J(ri,b  +  M,  V) 
=  b(t],o  -  Kn,^  -  O  + 

<\B{r],0\  +  \m,^-0\  +  \mr))\ 

if  the  interpolant  is  chosen  as  in  Lemma  4.1. 
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Moreover,  results  from  the  previous  theorem  provide  us  with: 

f 

I/K.X)!  <  C-^\\u\\,\mh,T,  <  C^I|a||,l|{||.,D,. 

<  C^\\uUV(\\o.T,  <  C^\\u\\,m\,.T„ 

SO  that 

IKII..®  £  c^llulk, 

and  this  completes  the  proof.  □ 

This  time,  the  rate  of  convergence  is  optimal  with  respect  to  h  but  the  rate  of  con¬ 
vergence  in  p  is  worse  than  in  the  previous  estimate.  This  makes  us  believe  that  the 
error  estimates  for  the  DG  method  can  still  be  improved  with  respect  to  p.  Maybe 
better  interpolants  are  yet  to  be  found. 


5.  Concluding  Remarks 

5.1.  Remarks  on  the  Discontinuous  Formulations 

We  have  studied  here  four  different  formulations  of  the  so-called  Discontinuous 
Galerkin  Method  (DGM).  These  formulations  simply  vary  by  one  sign  (plus  or  mi¬ 
nus)  and  by  the  addition  of  a  penalty  term  (or  not).  However,  they  greatly  differ  in 
nature  from  a  mathematical  point  of  view.  We  now  review  each  formulation  one 
by  one  and  recount  our  findings  in  the  case  of  linear  diffusion  problems. 

Global  Element  Method.  Little  can  be  proved  for  this  method.  We  were  able  to 
derive  tiie  continuity  of  the  associated  bilinear  form,  but  failed  to  even  obtain  a 
priori  error  estimates.  This  is  because  the  bilinear  form  is  not  guaranteed  to  be 
semi-positive  definite. 

Symmetric  Interior  Penalty  Galerkin  Method.  The  SIPG  Method  is  similar  to  the 
GEM  except  for  the  addition  of  the  penalty  term.  However,  it  allows  us  to  prove 
non  only  continuity  of  the  bilinear  form,  but  also  coercivity  in  the  discrete  discon¬ 
tinuous  space  (for  sufficiently  large  values  of  the  penalty  parameter),  and  thus  a 
priori  error  estimates  optimal  with  respect  to  h  (p,  —  1)  and  slighty  suboptimal  with 
respect  to  p  {s  —  3/2).  One  major  drawback  of  this  method  is  that  its  behavior  de¬ 
pends  on  the  selection  of  the  penalty  parameter.  If  not  chosen  carefully,  the  method 
can  fail. 
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Non-Symmetric  Interior  Penalty  Galerkin  Method.  The  limitation  of  the  SIPG 
method  is  remedied  by  changing  one  minus  sign  by  a  plus  sign.  Indeed,  although 
the  NIPG  formulation  results  in  a  non-symmetric  system  of  equations,  all  the  prop¬ 
erties  and  error  estimates  are  shown  to  be  independent  of  the  choice  of  the  penalty 
parameter.  We  also  find  the  same  rates  of  convergence  with  respect  to  h  and  p  as 
SIPG. 

Discontinuous  Galerkin  Method.  DGM  is  deduced  from  the  NIPG  method  by 
setting  the  penalty  parameter  to  zero.  We  then  observe  that  the  rate  of  convergence 
with  respect  tohorp  deteriorates.  Also,  in  the  case  of  the  pure  Laplacian  operator, 
when  c  is  set  to  zero  in  the  Poisson  problem,  we  obtain  a  priori  error  estimates  only 
by  defining  some  new  interpolants  whose  fluxes  are  weakly  equal  to  the  fluxes  of 
the  exact  solution  over  each  edge  of  the  elements.  Although  the  rate  of  convergence 
in  h  remains  optimal,  the  one  in  p  is  then  estimated  to  be  s  —  5/2.  We  believe  that 
it  might  be  possible  to  improve  this  rate  of  convergence  by  considering  other  types 
of  interpolants.  At  this  point,  detailed  numerical  experiments  would  be  helpful  to 
understand  how  the  penalty  term  affects  the  quality  of  the  approximations. 

5.2.  Future  Challenges 

The  great  challenges  for  IXlMs  are  to  1)  prove  uniqueness  of  the  solutions  of  the 
continuous  formulations,  2)  perform  more  numerical  experiments  to  understand 
the  role  played  by  the  penalty  terms,  3)  still  improve  the  a  priori  error  estimates  for 
the  Discontinuous  Galerkin  Method  of  Baumann  and  Oden,  4)  derive  rigorous  a 
posteriori  error  estimates  for  the  various  formulations. 


A.  Lemmas 

A.l.  Discrete  Schwarz  Inequality 

Lemma  A.l  Let  {a,}  and  {&,}  define  two  sequences  ofN  real  numbers.  Then 

N  /  N  \  /  N  \ 

X  ^  f  X  ( X  (A-i) 

Proof:  We  shall  show  the  discrete  Schwarz  inequality  for  N  =  2  first.  We  have: 

(flifci  -I-  a2b2f  =  a^bl  -I-  a^bl  +  2aibia2b2 

—  {al  +  «2)(^i  +  +  2aibia2b2 

=  (a\  -f  al){b\  4-  hi)  “  («i&2  -  afbif 
<  {a\  -I-  +  ^i) 
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so  that: 

H"  fl2^2  —  ^ “h 

The  result  is  easily  extended  to  N  >  2  by  recursivity.  □ 

A.2.  Multiplicative  Trace  Inequalities 

Lemma  A.2  Let  Q.  define  a  star-shaped  domain  with  boundary  dQ  as  shown  in  Fig.  4. 
Then,  for  all  V  e 

IWI§,.a<^( 


Cl 


■  sup|x| 
x^Cl 


|o||o,nl|Vz;||o,n 


(A.2) 


5Q 


P  |x|  <  x.n 


Figure  4.  Star-shaped  domain. 


Proof:  Let  O  €  Q  be  the  origin  and  let  n  denote  the  unit  normal  outward  vector  on 
dLl.  From  the  definition  of  a  star-shaped  domain,  there  exists  a  positive  constant  /3 
such  that 

;0|x|  <  X  •  n. 

Applying  Green's  Theorem  for  the  vector  field  u^x,  we  have: 

f  u\  •  nds  =  f  V  •  (m^x) dx. 

Jdci  Ja 

By  the  property  of  star-shaped  domains,  the  first  integral  is  shown  to  be  bounded 
below: 
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On  the  other  hand,  the  second  integral  is  bounded  above: 

f  V •  (u^x) dx=  f  -x  +  x-Vu^dx 
Ja  Jci 

=  [  2u^dx+  f  lux- Vudx 
Jci  Jci 

<2||M||o,i2+  f  \ux-Vu\dx 
JQ. 

2||w||o,Q  +  2sup|x|  f  |M||VM|dx 


< 


<  2||m||o,q  +  2sup  |x|||m||o,q||  VM||o,n 

X€il 


Using  both  bounds,  we  arrive  at: 

/ 

2 


ll«llo,dn  <  1 


inf  X 

xedii'  ' 


l|w|lo,Q  +  sup|x|||M||o,n||VM||o,Q 


which  completes  the  proof. 


□ 


Lemma  A.3  Let  K  be  a  triangle  or  a  quadrilateral  such  that  h^  <  gpj^  (shape  regular). 
Then,  for  all  veH^iK), 

Iloilo, 5K  ^  +  ||i;||o,k||  Vz;||o,Kj  •  (A.3) 

where  C  is  a  positive  constant. 


Proof:  Let  the  origin  O  be  the  center  of  the  inscribed  circle  in  K  with  radius  /Ok/2. 
We  therefore  have: 


sup  |x|  <  hK 

xGK 


W  N  >  pK>hK/Q 

x€dK 


SO  that  from  (A  .2) 


l“llo,aic^  j~  (|I“IIo,/<c  +  ^kII“IIo,k||Vm||o,k) 

^  2^  ^^||m||o,k  +  Iloilo, xIIVhIIo,/*:^ 


The  proof  is  complete  when  choosing  C  —  2q. 


□ 
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A.3.  Poincar^Friediich's  Inequalities 

Lemma  A.4  Let  Q  be  an  open,  bounded,  connected  domain  oj 
dCl.  Let  V  £  H^Q,  such  that 


Jo. 


vdx  =  0. 


Then 


lbllo,n  <  C||Vy||o,n 
where  C  =  C(Q)  is  a  positive  constant. 

Proof:  See  Schwab  [21,  p.350]  and  Brenner  and  Scott  [8,  p.l02]. 


with  Lipschitz  boundary 
(A.4) 

(A.5) 


□ 


Lemma  A.5  Let  z  G  Pp^(K)  and  z  be  the  average  ofz  onK,z  —  (/j^  zda:)/|X|.  Then 

i|z-z|kjc<C/zK||Vz||o,K  (A.6) 

where  C  is  a  positive  constant  independent  of  K  and  z. 

Proof:  Let  z  €  Pp^{K)  and  v  =  z  —  z.  Then 

f  vdx=  f  z  —  zdx=  f  zdx—  f  zdz  =  |fC|z  -  z|fC|  =  0. 

Ja  Ja  Jci  JQ 

By  a  scaling  argument  and  Lemma  A.4, 

Ikllo.K  <  CflK||t)||o,g  <  CWfeK||Vt;||o,K  <  ChK\\Vv\\o,K 

Substituting  z  —  z  for  v,  it  follows  that  ||z  — z||o,x  <  C/ik||V(z  —  z)||o,k/  in  other 
words,  since  z  is  constant,  ||z  —  z||o,x  <  C1ik||  Vz||o,x.  Q 

A.4.  Inverse  Property 
Lemma  A.6  Let  z  G  Pp^(fC).  Then 


IIVzllo,,  <  c£S||2||„,k 

«K 


(A.7) 


Proof:  See  Schwab  [21,  p.208]. 


□ 
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A.5,  Interpolation  Error  Estimates 

Lemma  A.7  Let  K  be  a  triangle  or  parallehgram  element  of  the  partition  2^  and  u  a 


^  -  - 

ofu,  Pk,  and  hf:,  and  a  sequence  Zp  €  Pp^,. 

(R),Pk  =  1,2,.. 

. ,  such  that  for  any  q,0<q  <s 

h'"~‘’ 

\\u  Zp||^,K<C 

Pk 

s  >  0 

(A.8) 

hi^-y2 

||m  -Zpilo,7  <C  ^  J  ,2  \\u\\s,K, 

Pk 

1 

^>2 

(A.9) 

where  p  =  min(pK  -M, s),  /ir  =  diam  (K)  and  7  C  dK. 

Proof:  See  Babuska  and  Suri  [3,4]. 

□ 

REFERENCES 

1.  D.  N.  Arnold.  An  interior  penalty  finite  element  method  with  discontinuous 
elements.  SIAM  J.  Numer.  Anal.,  19:742-,  1978. 

2.  I.  Babuska,  C.  E.  Baumann,  and  J.  T.  Oden.  A  discontinuous  hp  finite  element 
method  for  diffusion  problems:  1-D  analysis.  Computers  and  Mathematics  with 
Applications,  37:103-122, 1999. 

3.  I.  Babugka  and  M.  Suri.  The  h-p  version  of  the  finite  element  method  with 
quasiuniform  meshes.  Mathematical  Modeling  and  Numerical  Analysis,  21(2):199- 
238, 1987. 

4.  I.  Babuska  and  M.  Suri.  The  optimal  convergence  rate  of  the  p-version  of  the 
finite  element  method.  SIAM  J.  Numer.  Anal,  24(4),  1987. 

5.  G.  A.  Baker.  Finite  element  methods  for  elliptic  equations  using  nonconform¬ 
ing  elements.  Mathematics  of  Computations,  31(137):45-59, 1977. 

6.  G.  A.  Baker,  W.  N.  Jureidini,  and  O.  A.  Karakashian.  Piecewise  solenoidal  vec¬ 
tor  fields  and  the  Stokes  problem.  SIAM  J.  Numer.  Anal,  27(6):1466-1485, 1990. 

7.  C.  E.  Baumann.  An  h-p  adaptive  discontinuous  finite  element  method  for  computa¬ 
tional  fluid  dynamics.  PhD  thesis.  The  University  of  Texas  at  Austin,  1997. 

8.  S.  C.  Brenner  and  L.  R.  Scott.  The  Mathematical  Theory  of  Finite  Element  Methods. 
Springer-Verlag,  New-York,  1994. 

9.  Z.  Chen.  On  the  relationship  of  various  discontinuous  finite  element  meth¬ 
ods  for  second-order  elliptic  equations.  SMU  Math  Report  2000-02,  Southern 
Methodist  University,  2000. 

10.  B.  Cockbvurn,  G.  E.  Kamiadakis,  and  C.  W.  Shu,  editors.  Discontinuous  Galerkin 
Methods  -  Theory,  Computation  and  Applications,  volume  11  of  Lectures  Notes  in 
Computational  Science  and  Engineering.  Springer,  Berlin,  2000. 

11.  L.  M.  Delves  and  C.  A.  Hall.  An  implicit  matching  principle  for  global  element 
calculations.  /.  Inst.  Math.  Appl,  23(2) :223— 234, 1979. 


A  Priori  Error  Estimation  for  DGM 


41 


12.  L.  M.  Delves  and  C.  Phillips.  A  fast  implementation  of  the  global  element 
method.  /.  Inst.  Math.  Appl,  25(2):177-197, 1980. 

13.  J.  A.  Hendry  and  L.  M.  Delves.  The  global  element  method  applied  to  a  har¬ 
monic  mixed  botmdary  value  problem.  /.  Comput.  Phys.,  33:33-,  1979. 

14.  J.  A.  Hendry,  L.  M.  Delves,  and  C.  Phillips.  Numerical  experience  with  the 
global  element  method.  In  Mathematics  of  finite  elements  and  applications,  III, 
pages  341-348.  Academic  Press,  London-New  York,  1979.  (Proc.  Third  MAFE- 
LAP  Conf.,  Brunei  Univ.,  Uxbridge,  1978). 

15.  P.  Houston,  C.  Schwab,  and  E.  Suli.  Discontinuous  /ip-finite  element  methods 
for  advection-diffusion  problems.  Technical  Report  no.  00/ 15,  Oxford  Univer¬ 
sity  Computing  Laboratory,  2000. 

16.  O.  A.  Karakashian  and  W.  N.  Jureidini.  A  nonconforming  finite  element 
method  for  the  stationary  Navier-Stokes  equations.  SIAM  J.  Numer.  Anal, 
35(1):93-120,  1998. 

17.  J.  T.  Oden,  I.  Babuska,  and  C.  E.  Baumann.  A  discontinuous  hp  finite  element 
method  for  diffusion  problems.  TICAM  Report  97-21,  The  University  of  Texas 
at  Austin,  1997. 

18.  J.  T.  Oden,  I.  Babuska,  and  C.  E.  Baumann.  A  discontinuous  hp  finite  element 
method  for  diffusion  problems.  /.  Comput.  Phys.,  146:491-519, 1998. 

19.  B.  Riviere.  Discontinuous  Galerkin  Methods  for  Solving  the  Miscible  Displacement 
Problem  in  Porous  Media.  PhD  thesis.  The  University  of  Texas  at  Austin,  May 
2000. 

20.  B.  Riviere,  M.  F.  Wheeler,  and  V.  Girault.  Improved  energy  estimates  for  interior 
penalty,  constrained  and  discontinuous  Galerkin  methods  for  elliptic  problems. 
Part  I.  Computational  Geosciences,  3-4:337-360, 1999. 

21.  C.  Schwab,  p-  and  hp-  Finite  Element  Methods  -  Theory  and  Applications  in  Solid 
and  Fluid  Mechanics.  Oxford  University  Press,  New  York,  1998. 

22.  E.  Siili,  C.  Schwab,  and  P.  Houston.  /ip-DGFEM  for  partial  differential  equa¬ 
tions  witii  nonnegative  characteristic  form.  In  B.  Cockburn,  G.  E.  Kamiadakis, 
and  C.  W.  Shu,  editors.  Discontinuous  Galerkin  Methods,  volume  11  of  Lectures 
Notes  in  Computational  Science  and  Engineering,  pages  221-230.  Springer,  Berlin, 
2000. 

23.  M.  F.  Wheeler.  An  elliptic  collocation-finite  element  method  with  interior 
penalties.  SIAM  J.  Numer.  Anal,  15:152-,  1978. 


